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ABSTRACT 


Helicopter rotors operate in complex aerodynamic environment comprising of 
several significant fluid dynamic phenomena like blade-vortex interactions, reverse flow 
and stall on the retreating blade, formation of shockwaves and other compressibility 
effects on the advancing blade, and radial flow. This makes the rotor airloads calculation 
a difficult problem. Several aerodynamic models have been developed, starting from 
simple momentum theoiy, blade element theory, vortex theory to more complex CFD 
methods. Of these combined blade element-momentum theory is widely used for most of 
aeroelastic studies, because of its simplicity. However vortex methods are used for 
airloads computation. The aim of this study is to develop a prescribed wake vortex model 
to predict the performance of rotor blades. 

Earlier, a model based on the vortex theoiy was developed for rotor airloads 
calculation for a rectangular straight blade in hover. In the present analysis it is further 
improved to accommodate the effect of tip sweep. The effect of swept tip on 
aerodynamic performance characteristics are calculated and compared with those of 
straight rectangular blade results. It is observed that (aft sweep) reduction in rotor power, 
which has been observed in experimental studies. 



Chapter 1 
INTRODUCTION 


1.1. GENERAL 

A helicopter differs from an aeroplane in the sense that it can hover, take-off and 
land vertically, which makes it an extremely useful flying machine in times of distress or 
war. Because of the very demanding maneuvers it is expected to perform, the helicopter 
is a machine of considerable mechanical and aerodynamic complexity. This can be 
gauged from the slow progress in its development as compared to that of the aeroplane. 
Helicopters have come a long way from the noisy Gyroplanes and Sikorskys of the early 

tK 

20 century to the sophisticated Apaches. This could be made possible only through 
sustained research, which has led to improved performance, lift capability, and reliability. 
Since the 80s there has been an accelerated effort towards understanding and overcoming 
some of the technical, especially, aerodynamic problems in helicopter flight. With the 
rapid progress in processing speeds of computers, CFD techniques- mostly Euler and 
Navier-Stokes codes, are being applied to areas, which were, previously, unsolvable 
through traditional theories. 

Rotor airloads calculation is an aeroelastic problem. It involves the calculation of 
aerodynamic forces and the resulting rotor response iteratively. In forward flight, the 
difference in the velocities in the advancing and retreating sides cause a 1/rev variation in 
the aerodynamic environment, which means that the loading and the resulting response of 
the blades will be periodic. The 1/rev variation in velocities causes a rolling moment and 
hence very large 1/rev structural loads at the root of the blade. To counter these loads by 
aerodynamic forces, flap and lag motions are allowed. Pitch motion is allowed for blade 
control and elastic torsion motion occurs due to blade torsion flexibility. Elastic torsion 
motion changes the blade angle-of-attack, mad thus affects the airloads significantly. It is 



seen that, in general, the influence of structural dynamics on the calculated airloads is 
comparable to the influence of rotor wake effects, as is characteristic of aeroelastic 
problems. 

It can be seen that rotor airloads calculation is an extremely difficult problem and 
it is the belief of some researchers that not much progress has been made in this area. 
Bousman' delves into this aspect in a recent paper. He notes that, in rotor loads 
calculation, some of the present day analyses compute the magnitude of the critical 
design loads in maneuvers, correctly; but they fail to give the correct phase, which 
indicates that there is a deficiency in modeling of the physics involved. One such 
unsolved problem is the inability of the analyses to predict the negative lift in high-speed 
flight (Fig. 1). Another significant barrier to progress in the prediction of loads and 
vibration is the inability of the lifting-line analyses to predict the section pitching 
moments and hence the resulting control system loads (Fig. 2). 

The analytic method for rotor loads calculation can be summarized by the flow 
diagram given by Sadler^, Fig. 3. One set of the model deals with the wake model and 
associated calculations; the other deals with the determination of the aerodynamic blade 
loads and response. 

The flow field^ of helicopter rotors is extremely complex, as can be seen from 
Fig. 4. It comprises of several significant fluid dynamic phenomena like strong tip 
vortices, blade-vortex interactions, reverse flow and stall on the retreating blade, 
formation of shockwaves and other compressibility effects on the advancing blade, and 
radial flow. All these phenomena have significant effects on rotor performance, loads and 
vibration. For an accurate comprehensive analysis, all these features are to be included. 

Bousman*, in his paper, states that since the present day analytical methods do not 
consistently and accurately predict the loads and vibration of a rotorcraft, it is not 
possible to determine the relative importance of the various aspects involved in the 
modeling, like the rotor wake, unsteady aerodynamics, radial flow and elastic 
deformations of the blade. He, also, concludes that for the correct prediction of the loads 
and vibration, it is important that all the features be properly represented. 


1.2. ROTOR AERODYNAMICS 


The helicopter uses its rotating wings to provide for lift, propulsion and control 
forces. Major differences in the aerodynamic characteristics of a rotor from those of fixed 
wings are: 

1. High concentration of bound circulation over the outer portion of the blade 
resulting in an intense vortex trailed from the tip. 

2. A curved spiraling wake, which initially remains close to the rotor causing strong 
blade/vortex interactions. 

3. time dependency of stall and compressibility effects on the blade. 

Helicopter motion is of three types: 

1 . Hover : 

In this state, the rotor has no vertical or sideways velocity relative to 
air. 

2. Axial flight : 

In this state, the helicopter moves only vertically. 

3. Forward flight: 

In this state, there is, only, horizontal movement of the helicopter. 

Only hover state have been consider in the present analysis. 

1.2.1. Hover 

Figures shows the schematic of the wake for one blade, firom the flow 
visualization results of Gray‘d for hovering condition. As can be seen, the most dominant 
flow feature in the helicopter flow field are the strong vortices trailed at the tips of each 
blade which occur due to the concentrated lift and circulation there. Figured shows the 
visualization of the tip vortex shed^ by a blade. It is clearly showing the contraction and 
the convection of the tip vortex of one blade. 


The calculations of Clark and Leiper clearly show the wake feature® found 
experimentally by Landgrebe"*, namely, that the outer trailing vortices roll up quickly to 
form a strong tip vortex, while the inner vortices move downwards as a vortex sheet 
which becomes progressively more inclined to the rotor plane. Figure? shows the results 
of such calculations for the six-bladed CH-53A rotor. They also find that the initial 
position of the tip vortex depends strongly on the number of blades, if the number of 
blades is high, the tip vortex remains roughly in the plane of the rotor until it becomes 
close to the succeeding blade, when it is convected downwards. 

The other feature is the inboard vortex sheet of trailed vorticity, with opposite 
sign as the tip vortex. The axial velocity of the outer end of the inboard vortex sheet is 
greater than the tip vortex with the result that they separate. Also, the inclination of the 
inboard vortex sheet is due to the varying axial velocity along the radius. The aggregate 
of similar wake structures of all the blades of the rotor make up the complete wake 
representation. 

It is seen in Fig. 5 that successive turns of the tip-vortex remain a fixed distance 
away from each other. Also, the tip-vortex trajectory is taken to be a smoothly 
contracting semi-infinite cylindrical surface. This is, however, only a simplified 
illustration of the very complex nature of the rotor wake. Experiments have shown the far 
wake to be unstable"*. It was also seen that the tip- vortices from adjacent blades merge^’ 
Jain and Conlisk^ have shown this in their computations also. The interaction, which 
occurs in the region of maximum wake contraction, results in core breakdown followed 
by the diffusion of vorticity. This causes the wake to expand and hence there is a 
redistribution of vorticity in this region, which is close enough to the blade to induce 
significant inflow on the rotor. The said interaction was found to occur near the wake 
azimuth corresponding to the fourth tip vortex passage below the reference blade 
irrespective of the number of blades of the rotor. In order to make up for this additional 
source of inflow, Kocurek and Tangier’ proposed a vortex ring of radius equal to the 
rotor radius and at the axial level corresponding to the fourth tip vortex passage beneath 
the reference blade. Its strength was taken to be equal to four times that of the tip vortex. 
(Fig.8). In recent study the vortex model’® (Fig.9) used for calculation during pitch up 



motion. The vortex wake is represented by the group of concentric and coplanar vortex 
rings representing the evolution of the circulation along the blade with azimuth. 


1.3. INTRODUCTION TO TIP SHAPES 


The tips of the blades play a very important role in the aerodynamic perfomrance 
of the rotor. The blade tips encounter the highest dynamic pressure and highest Mach 
number, and strong trailed tip vortices are produced there. Figure 10 shows various blade 
tip designs” that have been used or proposed for helicopter rotors. There are several 
common designs, comprising those with taper, those with sweep, and those with a 
combination of sweep and taper. Some blade tips may also use an anhedral. Sweeping the 
leading edge of the blade reduces the Mach number normal to the leading edge of the 
blade, so allowing the rotor to attain a higher advance ratio before compressibility effects 
manifest as an increase in power required. However, the problem of rotor tip vortex 
formation and the effects of tip shape on the vortex characteristics such as velocity profile 
and diffusive characteristics is, however, still the subject of ongoing research. 

McVeigh and McHugh have conducted experiments with subscale rotors to 
study the effects of tip shapes on overall rotor performance and cruise lift-to-drag ratio. It 
was shown that the combined use of improved airfoil sections and tapered tip shapes can 
help minimize profile power and significantly improve rotor cruise efficiency. The effects 
of tip shapes were examined for four rotors having rectangular, swept, swept-tapered, and 
tapered tips. All rotors were tested at the same lift, propulsive force, and trim state, which 
provides a datum for performance comparisons. The tapered tip is found to give about 
10% higher equivalent L/D compared to the rectangular blade, but interestingly enough 
the rectangular blade gives a better maximum cruise LTD than either of the swept or 
swept tapered blades. This is because both the advancing and retreating blade 
characteristics of the tip shape are important, and an integrated performance metric such 
as L/D does not allow to distinguish separately between these characteristics. While 
sweepback alone clearly has the advantage of delaying the onset of compressibility 
effects to higher advanced ratios, the sweep may also promote early flow separation on 



the retreating blade at lower angle of attack. Therefore there can also be a performance 
penalty associated with a swept tip. 

Various types of tip shapes including rectangular, tapered, swept rectangular, and 
swept tapered were discussed thoroughly’^. Extensive research has been conducted in 
recent years to further quantify the benefit of various tip shapes. The following 
observations were found in the literature. 

Thomas H. Maier’"* studied the aeroelastic stability of straight and swept tip rotor 
plates. A planform view of the blades is shown in Fig .11. The swept tip blade has 30° of 
sweep which starts at 90% span. He found that for straight blades the flap, lag 
frequencies were slightly higher and torsion frequency was slightly lower than the swept- 
tip blades. 

The BERP (British Experimental Rotor Program) rotor blade is also a special type 
rotor because of its unique tip shape (Fig. 12.). This rotor was designed specifically to 
meet the conflicting aerodynamic requirements of the advancing and retreating blade 
conditions, either of which can limit the performance of the rotor in high-speed forward 
flight. However, the aerodynamic improvements shown with the BERP rotor are the 
result of several innovations in both airfoil design and tip shape design. One of the most 
recognizable features of the BERP blade is the use of high sweepback over the tip region. 
It reduces the compressibility effects and delays their effects on the rotor to a higher 
advance ratio. The BERP rotor is specially to perform as a swept tip at high Mach 
numbers and low angles of attack, but it is also designed to operate at very high angles of 
attack without stalling. The design of BERP rotor attempts to reduce compressibility 
effects on the advancing blade while simultaneously delaying the onset of retreating 
blade stall. 

Wind tunnel tests of a rotor were conducted’^ to determine the effects on dynamic 
response and aerodynamic performance by varying the tip design of the outboard 8% 
radius. Four different blade tip geometries or shapes having different planform, sweep, 
taper, and anhedral were tested. Results from the tests showed that blade torsional 
moments and control forces were reduced by adding sweep or anhedral. The anhedral tip 
benefited hover performance while the swept/taper tip provided the best performance at 
high advance ratios. 



Another wind tunnel test'^* was conducted to measure the aerodynamic benefit of 
a parabolic swept-back tip with anhedral on a 3-bladed model rotor equipped with very 
rigid blades. The swept onset was at 90% radius. Significant improvements on hover and 
forward flight performance and noise reduction over a rectangular planform were 
achieved. 

A test of a full scale rotor' ^ was conducted for the four tips (rectangular, tapered, 
swept rectangular, and swept/tapered) shown in Fig. 13, which were interchangeable over 
the outer 5% radius. All tips used a 9.5% airfoil. Data showed that the tapered planform 
had the largest impact on the overall performance while sweep effect was secondary. 

1.4. EXISTING THEORIES OF ANALYSIS 

The task of the aerodynamic analysis of the rotor blade is to find the induced 
velocity and performance characteristics for a given thrust. Aerodynamic models used in 
comprehensive analysis of the helicopter rotors have evolved from the simple momentum 
theory based on a uniformly loaded actuator disk through the classical blade element- 
momentum method and wake modeling procedures to CFD methods which developed in 
pace with the high speed digital computer. 

Momentum theory applies the basic conservation laws of fluid mechanics 
(conservation of mass, momentum and energy) to the rotor and the flow as a whole to 
estimate the rotor performance. It is a global analysis, relating the overall mean flow 
velocities and the total rotor thrust and power. The rotor is assumed to be composed of an 
infinite number of blades and operates in an incompressible fnctionless fluid. Momentum 
theory provides a fast, convenient means of predicting overall performance, but is not 
concerned with the details of rotor loads or flow and hence alone is not sufficient for 
designing the rotor blades. 

Blade element theory calculates the forces on the blades due to its motion through 
the air, and hence the forces and performance of the entire rotor. Basically, blade element 
theory is lifting-line theory applied to the rotating wing. Each section of the blade is 
assumed to act as a two-dimensional airfoil to produce aerodynamic forces, with the 
influence of the wake and the rest of the rotor contained entirely in the induced angle of 


attack at the section. The solution thus requires an estimate of the wake-induced velocity 
at the rotor disk, which is provided by momentum theory’, vortex theory, or combined 
blade-element/momentum theory. 

Vortex theory is a rotor analysis that calculates the flow field of the rotor wake 
and the induced velocity at the rotor disk by using the fluid dynamic laws governing the 
action and influence of vorticity (the Biot-Savart law, Kelvin’s theorem, and Helmholtz’s 
laws). This theory uses lifting line or lifting surface theories together with a wake model 
for the performance calculations. The wake geometry is taken to be either the reliable, 
easy-to-use prescribed wake or the sophisticated and computationally demanding 
wake. 

In the prescribed-wake analysis, the geometry of the vortex sheets from the 
individual blades are prescribed in advance, which implies that the velocity field has been 
assumed. In this approach, the wake geometry is specified as a frmction of rotor 
configuration and thrust level through simple analytical expressions. 

Landgrebe'^ demonstrated the practicality of this method by incorporating an 
experimentally derived generalized wake description in the UTRC Prescribed Wake 
Hover Performance Analysis. Kocurek and Tangier^ used a prescribed wake lifting 
surface model in their performance analysis. 

In the free-wake analysis, an initial distribution of the vortex sheets is assumed 
and the elements of the vortex sheets are allowed to convect in the velocity field they 
create. The vortex elements will move until they take up positions, which are consistent 
with the velocity field. This analysis requires orders of magnitude more computation than 
the prescribed wake. Application of the free-wake analysis can be found in the works of 
Sadler^, and Bagai and Leishman'®. 

Prescribed wake does not give good results in low speed forward flight while free 
wake does. But the importance of the free wake declines as airspeed increases and at hi^ 
speeds, there is no difference in the two wake models.’ 

Lifting line theory is based on the assumption that the wing has a high aspect 
ratio, which is, in general, valid for a helicopter rotor wing. However, in areas where the 
loading or induced velocity has high gradients such as the sections near the tip or near an 
encounter with a vortex from a preceding blade, lifting line theory may not be accurate. 


This model usually gives high gradients in induced velocity when the tip vortex passes 
close to the tip. Lifting line formulation neglects all planform effects. 

In second-order lifting line theory each blade is divided into a finite number 
of span-wise segments with a span-wise bound vortex located at the quarter chord line 
and a control point placed at the middle of three-quarter chord line in each segment. The 
wake vortex filaments are trailed from the trailing edge of the blade. Hence, there are the 
so-called chord-wise bound vortices between the lifting line and the trailing edge. In this 
assumption, the close encounter of the trailing vortices with the control points can be 
effectively avoided during the solution. Second order lifting line theory has given 
results same as lifting surface theory and can be used to improve the calculation of 
airloads instead of resorting to the computationally demanding lifting surface or CFD 
methods. 

Lifting surface analysis is required for an accurate treatment of the vortex-induced 
loads on a rotor blade, especially for extremely small vortex-blade separations where the 
loading varies greatly in a short distance along the blade. This theory represents the wing 
by vortex surfaces and satisfies the boundary conditions over the entire surface. It can 
handle the large variations of induced velocity and loading that occur at the blade tip or in 
an encounter with a wake vortex. In Kocurek and Tangier’s modef, the rotor lifting 
surface is developed from a “Vortex-Box “ surface and wake numerical model. 

Lifting surface analysis would be beneficial for any planform which varies 
significantly from the conventional, constant chord, blade. A sensible compromise 
between accuracy and computational effort would be to provide a single panel lifting 
surface model in rotor applications. The full lifting surface model may be used where a 
more detailed accounting of both planform and wake influences on the blade is needed. 

Baskin et al^^ give a very thorough exposition of the application of the vortex 
theory to rotary-wing aerodynamics in their book. Johnson provides a detailed 
comprehensive analysis of rotor wings based on the vortex methods and also an overview 
of some of the methods used for a comprehensive analysis. 

Potential theory finds the induced velocity by solving the fluid dynamic equations 
for the velocity or stream function. It requires the assumptions that the flowfield is 
inviscid, irrotational, and isentropic. Application of velocity and acceleration potentials 



makes it possible to determine steady and unsteady flow fields induced by the rotor in 
both incompressible and compressible fluids with a precision similar to that offered by 
the vortex theory approach but with less computational effort. Mangier and Squire were 
one of the first to adapt the velocity and acceleration potential concepts to the rotor for 
finding the induced velocity field. Rotary wing Aerodynamics^^ gives a brief outline of 
the application of potential theoiy to the determination of flow fields around three- 
dimensional, non-rotating bodies. 

Some of the current methods of used for calculating rotor performance 

solve the potential, Euler or Navier-Stokes equations coupled with an external free- or 
prescribed-wake model based on the lifting line or lifting surface theory. These 
approaches require fairly large computer resources from solving two coupled models 
simultaneously. The other approaches are the CFD wake capturing methods, wherein the 
entire flow field is simulated with a unified differential method. No wake model is used, 
provided the grid contains the whole domain. For the rotor in hover, CFD methods have 
been used for Euler computations or Navier-Stokes computations with different strategies 
like embedded grids or unstructured grids. In principle, CFD methods may be used to 
compute fundamental fluid properties on both the rotor blades and the wake. However, 
the degree of complexity associated with the CFD schemes demand prohibitively high 
computational tools. 

1.5. OBJECTIVES OF THE PRESENT STUDY 

The objective of this study is to improve the vortex model (having prescribed 
wake structure) developed in Ref. [29] for rectangular blades, to include the effect of tip 
sweep. In order to validate the mathematical formulation and the computer code, first the 
aerodynamic loads on a swept fixed wing were calculated. These results are compared 
with the results available in the literature. Then performance of rotor blades having swept 
tip have been evaluated and the effects of tip sweep on rotor thrust and rotor torque have 
been analysed. 


Chapter 2 

MATHEMATICAL MODELLING 


2.1. DESCRIPTION OF THE PROBLEM AND THE METHOD CHOSEN 

The aim is to develop a code to carry out the performance calculations of a rotor 
blade with swept tip in hover state, given the rotor configurations like the number of 
blades, size of the blades, pitch, sweep. 

Lifting line theory is used for the analysis considering the conclusion of the 
authors : ‘...the lifting line approach is adequate for predicting the hover performance of 
a wide range of conventional and advanced rotor designs.’ 

The wake geometry is specified based on the semi-empirical model hypothesized 
by Gray and further developed by Landgrebe"^. The equations of Landgrebe have been 
used as such but for the tip vortex parameters, which have been modified by Kocurek and 
Tangier^, who found Landgrebe’s generalized equations to be in disagreement with the 
tip vortex trajectories in their wake geometry investigations. Hence, in the final form, the 
tip vortex coordinates are given by Kocurek and Tangier’s recirculation model, and the 
geometry of the inboard wake sheet is given by Landgrebe’s wake model. The wake 
geometries are functions of the rotor thrust coefficient, number of blades, collective pitch 
angle, and blade twist. 

A rotating coordinate system, fixed to the rotor and moving with it is chosen 
(Fig. 14). Since for hover, the flow field is axisymmetric, it would be sufficient to find the 
loads on one blade. The results hold for all the other blades. The blade section is 
considered as 2-D infinite section. Conventional strip theory based on 2-dimensional 
airfoil data is assumed applicable to compute the rotor performance characteristics. 
Airfoil data has been taken from standard tables^^ for a A64C4 0012 blade section. 
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2.2. ASSUIVIPTIONS 


Assumptions made in the analysis are: 

1 . A rigid rotor is assumed with constant chord case was considered. 

2. Blade dynamics is not included. 

3. The rotor blades are given only small pitch angles. 

4. The wake is assumed to be made of line vortices, which have no core. 

5. Stall and compressibility effects have not been considered. 

6. Fluid dynamic characteristics of the vorticity like vortex-vortex interactions; 
stability of vortices; vortex core bursting and other viscous effects have not been 
taken into account. 

7. The trailing vortices from the swept tip are also taken parallel to the trailing 
vortices from the straight portion of the blade. 
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2.3. FORMULATION 


2.3.1. Coordinate Equations Used 

Tip Vortex Coordinates (Fig. 15a) 

The generalized empirical equations^ (Kocurek and Tangier’s model) defining the 
tip vortex axial and radial coordinates as a function of azimuth are given by. 

Axial Coordinates 

'^71 

2 , = K 0 < Xj/^. <xi/^ ,y/^ =— 

( 1 ) 

where Z/ is the axial coordinate of the vortex 


Radial Coordinates 

r = A+().-A)e~^'*''^ 0< ip <4y/ (2) 

where A is approximately 0.78. 

ki rate at which the vortex settles before it encounters the following blade 
kj rate at which the vortex settles after it passes the following blade 
X wake contraction parameter 

ki and kt are shown in Fig. 15a 


The tip vortex trajectory equations represent the wake’s structure upto the point of 

maximum contraction, which occurs close to a wake azimuth of — . 

N 

Vortex Settling Rate Parameter, h and k 2 



(3a) 



where 


B = -0.000729^, 

C = -2.3 + 0.206^, 
m = 1.0-0.25^°'”°''' 
n = 0.5-0.0172^, 

0^ total twist in the blade (washout) in degrees 

Ct Thrust Coefficient 

N number of blades 

k,=-iCr-C,J^ (3b) 

where 

Cr^ = N^i-B/Cy 


Wake Contraction Parameter, X 

I 

X = 4.0(Cry (4) 


Features; 

• the axial velocity of tip vortex element is very low until it passes beneath the 
following blade. At that point, the tip vortex element lies radially inboard of the 
following blade and the induced effect of the bound and trailing vortices of that 
blade increase the axial displacement rate. Hence, k 2 is greater than kj. 

• axial displacement and radial contraction of the tip vortex increase with increase 
in thrust level. 

• For twisted rotors, the tip vortex passes above the rotor disk at low values of Ct 
since the tip region is producing thrust in opposition to that produced further 
inboard. As Cj increases to where k^ = 0, there is no tip vortex and hence the tip 
becomes imloaded. The vortex then reappears as Ct is increased further. 

• k 2 also goes to zero when kj =0 



Vortex Sheet Coordinates (Fig. 15b) 

The equations for the inboard vortex sheet’s axial and radial coordinates have 
been used witliout alteration from Landgrebe’s model.'^ 

Axial Coordinates 

The inboard vortex sheet extends from the root cutout to the point where the tip 
vortex emerges. The axial coordinates of points on the cross sections of the vortex sheet 
shed by the inboard portions of the blades were approximated as varying linearly with the 
radius. These linear approximations are extended at both ends until they intercept the axis 
of rotation (F=0) or an imaginary cylinder of radius (F=l). Landgrebe gives the 
coordinates of these two end points. (Fig. 15b). 
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where the vortex sheet parameters are : 
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Radial coordinates : 


r = ( 6 ) 

where 

normalized radial location of the point on the blade where the 

T 

vortex originates = 

fj normalized radial location of the tip vortex where it has the same 
vertical displacement as that of a point on the vortex originating 
from 


2.3.2. Biot-Savart law: 


The Biot-Savart law is the basic relationship used for calculating induced flow in 
the vortex theory. For a line vortex with circulation F, the differential velocity induced, 
dv, by a differential length of the vortex, dl, is given by: (Fig. 16) 


r dJ xr.^ 

~^~W 


( 7 ) 


Sign Convention 

Looking in the direction of the vortex, if the rotation of vorticity is 
clockwise, then F > 0; if the rotation is anti-clockwise, then F < 0. 



2.3.3. Derivation for tip-vortex induced velocity expressions 

Referring to Fig. 17, 

Parametric equations of the tip-vortex helix 
=-rsin^„ 

y^=rcosip,^. ( 8 ) 

where 

r is defined in eq. (2) 

y/^. is measured clockwise since blade rotation is anti-clockwise 

Parametric equation of the lifting line of the reference blade 
X = - /? sin 0 

y = pcosQ (9) 

z = 0 

where p is the radial location of the reference point on the blade. 

Note: 

The angles have been taken zero because the lifting line of the reference 
blade is the line from which angles are measured. 

r,=P-R, (10) 

where Vectorial position of a point on the tip vortex 
P Vectorial position of a point on the lifting line 

Fj Vectorial position of a point on the lifting line with respect to a point 
on the tip vortex. 

P ■=x i +y j+z k 
where x,y,z given in eq. (9) 


where Xtv ,ytv, Ztv given in eq. (8) 

-> r, =(-/7sin0-.r„.)f+(/9COs0-:v„)j+(-z„,)^ 


(-psin 0 - ) ' + (pcosO - ) ' + (-z^ ) ' 


Applying Biot-Savart Law 


r, X H 
dv = — ^ 


An 


where dv is the velocity due to the infinitesimal element at any point, T, on 
the tip vortex 

the tangential vector at T, dJ , is given by 

if 7 

dl =- — .d\f/^. = 
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2.3.4. Derivation for Inboard Vortex sheet induced velocity expressions 


The inboard vortex sheet originates from points on the blade lying between the 
root cutout (F = 0.1) and the point at which the tip-vortex starts (F = 0.97 ). The cross 
section of the vortex sheet at each azimuth is taken to be linear.(Fig.8). Equation (5) gives 
the axial displacement coordinates of the extended vortex sheet, at the points F = 0 (axis 
of rotation) and F = 1 . From eq. (5), the axial coordinates for the intermediate points on 
the vortex sheet coordinates can be interpolated. 
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2.3.5. Derivation for far wake induced velocity expressions 

The far wake is represented by a vortex ring at a vertical distance corresponding 
to the fourth tip vortex passage below the reference blade. It is shovra in Fig. 18. 

From Biot-Savart law, 

j- Ty J/ X r 

dv=— 

An Iz-j^ 

where 

Fj- is the circulation of the far wake 

dv is the velocity induced at P due to the 

infinitesimal element at a point, T, on the vortex ring. 



where 


r =P-R, 

r Vectorial position of a point on the lifting line with respect to a 
point on the vortex ring 

P Vectorial position of a point on the lifting line 
Vectorial position of a point on the vortex ring 


.R, =-Rj- sin^i +Rj- cos^ j + hk 

where 

Rf radius of the far wake ring 

h vertical distance between the blade and the far wake 

— A /N 

P = - yOsinOz + pcoz 0 j+Qk 

— dR, A ^ 

dl =-^d^ = (-Rj-cos^i - R^sin^j+0k)d^ 

where (j) azimuthal angle of the far wake vortex element with respect to the 
reference blade 


dl xr = 


i 

-Rj-cos^ 
(-Rf sin^ + psinO) 


J 

-Rf sin^ 

(Rf cos^-pcosO) 


A 

0 

-h 


d(f) 



Of. 



2.3,6. Derivation for the induced velocity expressions due to the other lifting lines 

The bound vortices of the other blades too contribute to the induced velocity of 
the reference blade. (Fig. 19) 

From Biot-Savart Law, 

r(/) dixr 

where dl vectorial element of the other bound vortex 
r = P-P, 

where r is the vectorial location of a point on the reference lifting line with 
respect to an element of the other bound vortex 
P Vectorial position of a point on the lifting line 

P, Vectorial location of a point on the bound vortex 

P = - ps'mOi+pcosOj+Ok 

P, =~ Pi sm(0 + k* ST/2)i +pi cos(0 + k * / 2)J + Ok 

where k blade number 0,1,2....N; 0 being the reference blade 
Elemental length, dl, of the bound vortex is given by 

I = 0 z + Pi _/ + 0 P 
dl = 01 +dpij+0k 

A A 

i j k 

dIxr = 0 1 0 dpi 

(-/?sinO + p, sin(0 + A:*;Y’/2) (yOCOsO-p, cos(0 + it*;7r/2) 0 

V 

where rc is the root cutout and R is the radius of the blade 


r(/) dlxr 

Atc \r? 


2.3.7. Lifting Line Theory expressions 


From Fig. 20, it is seen that 


or, = G, -P, 

where 


00 

A. 

subscript j 


effective angle of attack 

total angle at a section 

Pitch of the blade 

twist at the section 

induced angle at the section 

number of the control point at the section 


( 12 ) 


From Joukowsky’s theory, 

l = pWjr. (13) 

where 

/ lift per unit length 
Wj resultant air velocity at the section 
Fj circulation at the section 

But, fi-om experimental aerodynamics, 

l = ^c,pw;c, (14) 

where 

Qj coefficient of lift at the section 
p density of air 
Cj chord at the section 




where 

Vy vertical component of the total induced velocity 
Vxj horizontal component of the total induced velocity 
perpendicular to the blade 

Uj relative velocity of air due to rotation of the blade = Q, rj 
Voo axial flow velocity 

Eq (16) is the final form of the lifting line equation used. This is solved for /) in 
an iterative process. 




2.3.8.Performance Parameters Derivation (Fig. 20) 


Consider an elemental blade segment A Fj . As shown in Fig. 16, 
The tlrrust per unit length. 


But, 


F,j = Lj cos -Dj sin /3j 
-F.j =L^sinfij+DjCos/3j 

Lj=C,.^Wj%Ar. 


( 17 ) 


(18) 


D- =Cr,-^W^b.Ar 

j 2 j j 

where Qy is the lijft coefficient and Cdj is the drag coefficient of the airfoil 


It can also be seen in Fig.20 that 
W- cos Pj = Hr J- v^j 

(19) 

IF,sin/?,. =F + v,, 

J * J y 


substituting (15), (18), and (19) into (16) and taking// ■ = 


Fy = Vj - v^j ) - (F + v^j )]rjAr. 


- Fy = PW + ) + Pj 


Torque generated by the segment 
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DJ 

^ c 

'-'iy 


AMj =-F,j.rj 


= PW + + Ay (-^y - v^.OFyO'^'y 



Total Thrust 


T = NZF^. 
Total Torque 


M = NZAM. 


Coefficient of Thrust, 
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Chapter 3 

SOLUTION METHODOLOGY 


3.1. METHODOLOGY 

The wake analysis involves two tasks - calculation of the wake-induced velocities 
and calculation of the resulting aerodynamic loading on the rotor blade. 

Basically, the wake-induced velocity on the rotor blade is obtained by integrating 
the Biot-Savart law applied to the vortex elements in the rotor wake. For convenience in 
calculations, the wake is divided into several components - tip-vortex, inboard vortex 
sheet, and the far wake, which is represented by a ring vortex. There is also the 
contribution of the lifting lines of the blades other than the reference blade. The extent of 

Stt 

the wake is taken to be upto an azimuth angle of — . The wake is divided into azimuth 

step-sizes of 15°. The integration is done using an adaptive Simpson’s rule (Appendix B) 
over these 15° steps. 

To obtain the blade loading, the blade bound circulation variation is to be 
determined. This is given by the solution of Prandtl’s lifting line theory. However, 
Prandtl’s lifting line equation, being an integro-differential equation^^ is difficult to solve. 
The lifting line equation is converted into a set of algebraic equations by applying the 
lifting line equation to only a prescribed set of control points. 

The wake from each blade is represented by a finite number of discrete vortex 
filaments - a strong, rolled up tip vortex filament, and several weaker filaments 
representing the inboard portion of the vortex sheet. The set of control points, mentioned 
above, are taken to be lying in between the points of origin of the vortex filaments (Fig. 
21). Since the blade circulation has a larger gradient near the tip, the control points in this 
region are placed at closer distances than in the rest of the blade. The blade bound 
circulation and induced velocities are calculated at these control points. The positioning 
of the control points in between the vortex filaments also ensure that singularities do not 
occur (Fig.21). The bound circulations at the control points are the unknowns in the set 



of algebraic equations obtained by applying the lifting line equation at the control points. 
These are solved by an iterative procedure. The difference between the bound 
circulations at the two adjacent control points defines the strength of trailed vorticity, 
which emanates between these points. Since the unknown circulations are solved 
iteratively, an initial estimate of the circulation would be required. For this, combined 
blade-element/momentum theory is used. The drop in circulation at the tip is, however not 
predicted by this theory, which is, hence, accounted for by setting the circulation to zero 
over a small distance at the tip. 

To compute the rotor performance characteristics, 2-dimensional airfoil data is 
assumed to be applicable. Lifting line theory does not predict the vortex roll-up process. 
However, it has been approximately modeled^^ by letting the discrete tip vortices to 
combine at an azimuth angle of 60°. The vortex wake constructed using empirical 
equations (section 2.3.1). The wake diagram shown in Fig.23. 



3.2. ALGORITHM OF THE CODE 


The algorithm for the code is given here and the flowchart of the program 
is given in Fig.22. 

1 . Rotor data is input. 

2. An initial value of the blade bound vortex is estimated using the combined blade- 
element/momentum theory or from the value of blade circulation calculated 
previously for rotor configuration close to the present configuration. The strength 
of the trailing vortices is calculated from the blade bound circulation. 

3. Next, the velocity induced at the blade by each of the wake component - tip 
vortex, inboard sheet, other bound vortices and the far wake- is calculated by 
applying the Biot-Savart law. 

4. The effective angle-of-attack at each section of the blade is then calculated. Using 
eq. 16 (section 2.3.7), the new values of blade circulation are obtained. 

5. If the new value of blade circulation is not comparable to the old value, blade 
circulation value is slightly increased or decreased, depending on the new value of 
the blade circulation. 

Note: in the present analysis, the new and old values of the blade circulation are 
taken to be comparable if their difference is less than 0.1 times the old 
value of circulation. Since it is taking more time for each iteration, this 
assumption was made. 

6. Steps 3-5 are iterated till convergence of blade circulation values. 

7. Now, firom the converged values of blade circulation and induced velocity values, 
the effective angle-of-attack, coefficient of thrust and coefficient of torque are 
calculated. 



Chapter 4 

RESULTS AND DISCUSSION 

4.1. VALIDATION FOR SWEPT FIXED WING 

The code was, initially, verified for the simpler case of swept tapered fixed wing 
since the wake consists of semi-infinite straight trailing vortices. The code was slightly 
modified for the swept fixed wing case but the base algorithm remained the same. The 
wing details are as follows; 

Wing airfoil section : NACA 65-210 
Aspect ratio, AR : 9.02 
Taper ratio : 0.4 

The code results were validated with lifting surface results^® for various sweep 
angles (as shown in Fig. 25) for different angle of attack. While increasing the angle of 
attack from 5“ to 10°, the results were found close to lifting surface results. The lift curve 
for the swept tapered fixed wing was obtained and compared with the straight fixed wing 
results. As shown in Fig. 26, the experimental data^° and the lift curve obtained from the 
code match to large extent, except near stall and it is also observed thast the lift for swept 
tapered wing was slightly decreasing as compared to straight wing. 

4.2. CALCULATIONS FOR THE RECTANGULAR ROTOR BLADE 

The vortex model^^ (having prescribed wake structure) was developed for 
rectangular blades. The following parameters were studied. 

1. Root cut-out; The lifting portion of the blade generally starts at a radial station of 
about 10% to 30% of the blade radius. The root cut-out was assumed as 10%. 
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2. Tip vortex factor (tvf): The extent of the tip region of the blade from which the tip 
vortices emerge (denoted as tip vortex factor, tvf). To fix the tvf, calculations for the 
different cases of the tvf were carried out and the value of 97% was proposed as it 
gave results closest to experimental data. 

3. Roll-up location: After several attempts, the roll-up of the tip vortices has been chosen 
to occur at 60° azimuth. 

4. Strength of the vortex ring: The strength of the vortex ring, which is added to make up 
for the additional inflow due to the diffused vorticity in the region below and 
outboard of the well-defined near wake. The vortex ring strength was proposed 8 
times that of the tip vortex strength, for which the results are in good agreement with 
the experimental data. 

The values proposed for various parameters 

1 . Root cut-out 10% 

2. Tip vortex factor (tvf) is 97 % 

3. Roll-up location of tip vortices at 60° 

4. Strength of vortex ring is 8 times strength of tip vortex 

were used in present study. Using the above parameters, the performance calculations 
were carried out for a NACA 0012 rotor wing (Fig 24). The other details of the model 
blade configuration^ are given in tablel. These rectangular blade results were validated 
with experimental data, as shown in Fig 27. It is observed that for small pitch angles the 
results were matching with experimental data but for large pitch angles (>10°) this niodel 
not predicting good results. 


UNTWISTED ROTOR 


Aspect Ratio 

N 

R 

m(in.) 

b 

m (in.) 

OR 

m/s (j^s) 

Rotor 

solidity 

(tj) 

18.2 

4 


.037(1.47) 

152.4(500) 

0.070 


Table 1. Rotor Configuration used for calculations (Ref. 7) 







4.3 CALCULATIONS FOR SWEPT TIP ROTOR BLADE 


The sweep angle of the blade starts from 95% rotor radius. Performance 
calculations were carried out for swept tip rotor blade. These results are compared with 
straight rectangular blade results, to identify the effect of tip sweep. 

Keeping the pitch angle 6”, the coefficient of thrust and torque were calculated for 
sweep angle varies 0° to 1 5°. The results are shown in table 2. (See Fig. 28). It is 
observed that the reduction in torque is more than reduction in thrust. 


N = 4, pitch angle = 6°. 


Sweep angle 

Cq / cr 

Ct/a 

0° 

0.00290 

0.0349 

T 

0.00235 

0.0348 

5° 

0.00230 

i 

0.0342 

15° 

0.00220 

0.0336 


Table 2. Effect of tip sweep on Torque and tluust. 

Keeping the sweep angle 5“, the performance calculations were carried out for 
various pitch angles. These results were compared with straight rectangular blade results, 
as shown in Table 3. Results indicate that there is more reduction in torque in comparison 
to thrust. 


Pitch angle 

Cq/a 

J 

Ct/ a 

Straight blade 

5° tip sweep 

Straight blade 

5° tip sweep 

4° 

0.0018 

0.0016 

0.019 

0.0187 

6° 

0.0029 

0.0023 

0.049 

0.0342 

8° 

0.0046 

0.0036 

0.053 

0.0521 

10° 

0.0067 

I 

0.0042 

0.0731 

0.0724 

12° 

0.0092 

0.0068 

0.0942 

0.0932 


Table 3. Comparison of swept-tip results with strai^t blade results. 
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Figure 29 shows the effect of 5° tip sweep of rotor blade on performance. It is observed 
that (aft sweep) reduction in rotor power, which has been observed in experimental 
studies. 

Keeping the pitch angle 6°, the effects of tip sweep on bound circulation, 
effective angle of attack, induced velocity, induced angle of attack and load distribution 
along the span were calculated. The effect of tip sweep on blade boun circulation is 
shown in Fig 30. The reduction in bound circulation was observed. Figure 3 1 shows the 
effect of tip sweep on the effective angle of attack. The results showed that the reduction 
in angle of attack as compared to rectangular blade results. The load distribution along 
the span shown in Fig 32 for various sweeps angles. The reduction in lift was found as 
similar observation was made for fixed swept wing. The variation of induced velocity 
along the span for various tip sweep angles was shown in Fig 33. The increse in induced 
velocity was found as comparatively rectangular blade. The induced angle of attack along 
the span for various tip sweep angles was shown in Fig 34. The results shows, near the 
blade root the increse in induced angle of attack is more as comparatively tip. 



Chapter 5 
CONCLUSIONS 

In this study, a prescribed wake model has been developed to study the hover 
performance of a swept tip rotor blade. The model has been used to evaluate the rotor 
thrust and rotor torque for various values of tip sweep and rotor pitch angles. The 
important observation of this study is that tip sweep reduces both rotor thrust and rotor 
torque in comparison to a straight blade. The reduction in torque is observed to be more 
pronounced than the reduction in thrust. A similar observation has also been made in 
experimental rotor tests. 
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Chapter 6 

SUGGESTIONS FOR FUTURE WORK 

1 . In this study, the rectangular swept tip blade was considered. The model can be 
improved by including taper effect also. 

2. Additional aerodynamic phenomena related to compressibility, radial flow and 
dynamic stall effects may be included to develop a prescribed wake model 
applicable for forward flight. 



APPENDIX A 


Combined Momentum /Blade element Equations 
{used for non-uniform inflow calculation) (Ref 31) 



dCr =^{dr^-Zr)dr 


where 

axial velocity 
Vi induced velocity 
Q angular velocity of rotor 
G solidity ratio 

6 total angle of attack at the section 
a lift curve slope 
R radius of the rotor blade 
Ct coefficient of thrust 


APPENDIX B 


Simpson ’s Rule of Integration (Ref. 32) 

Algorithm 

b 

This algorithm computes the integral J = ^ /(^') dx from given values 

a 

fj = fi^j) equidistant xo = a , xi = xo + h , X 2 m = xq + 2mh = b by 

Simpson’s rule where h = (b-a)/2m 

INPUT : a,b,m,fo ,f 2 m 

OUTPUT : approximate value J of J 

Compute So = fo+f 2 m 

51 = fl+f3+ +f2m-l 

52 = f2+f4+- • • •+f2m-2 

h = (b-a)/2m 
J = hy3(So+4Si+2S2) 

OUTPUT J Stop 
End SIMPSON 

Adaptive Integration 

Here, the step h is adapted to the variability of f(x) : small h, where the 
variability of f(x) is large and large h where it is small. Changing h is done 
systematically, usually by halving h, depending on the size of the estimated error over a 
subinterval. The subinterval is halved if the corresponding error is still too large, that is, 
larger than a given tolerance(maximum admissible absolute error) or is not halved if the 


error is less than the tolerance. 
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Fig. 1 , Airloads Estimations for Two Radial Stations (Ref. I) 
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Blade loads and 
response program 


Fig. 3. Flowchart of Sadler’s Model (Ref. 2) 
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Fig. 5. Gray’s Schematic of Hovering Wake Structure (Ref. 4) 



Fig. 6.Visualization of Tip Vortex Shed (Ref.5) 
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Fig. 7. Vortex Wake Calculated from Free-Wake Analysis (Ref.6) 



Fig. 8. Cross-Section of Wake With Vortex-Ring (Ref.7) 
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Fig.lO.Blade Tip Shapes (Ref.l 1) 





Figure. 12.BERP Rotor Blade (Ref. 11) 
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Fig. 13. Tip planfonn geometry (Ref. 13) 
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axial coordinate, - ij 



Fig. 14. Rotor Coordinate System. 
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Fig. 15a. Tip Vortex Trajectory. (Ref.5). 
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Fig 15b. Inboard Vortex Sheet Coordinates. (Ref.4). 
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Fig. 18 . diagram showing the position of the vortex ring. 
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Fig,20. Flow at a Blade Section. 
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Fig. 21. Trailing Vortices From the Swept-Tip Rotor Blade 
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Fig. 22 . Flow Chart of the Program 
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Fig. 22. Flow Chart of the Program 
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Woko iTom cxYs ot tt>o Wadss crt a twob^aded rotor in hover 
(pitch » 4 deg.) 



Fig 23a. Wake Constructed Using Empirical Equations 


Sicle-Vi€w 



Fig 23b. Side view of the Wake 
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Fig 24. NACA 0012 Curves 
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Coefficient of Lift, CL CL/( Aspect ratio) per degiec 



Fig. 25 Finte wing validation (The numbers by the symbols represent 
sweep angles in degrees) 



Fig 26. Effect of sweep on fixed wing lift 
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Effective Angle of Attack(degrees) ^ Blade ^und circulation 


Effect of Tip Sweep on Circulation 



Effect of Tip Sweep on Angle of Attack 
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Load per unit length (kg/s) 



r/R (normalized radius) 

Fig 32. Variation of Load along the span 
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Induced velocity (m/s) 



r/R (normalized radius) 

Fig 33, Variation of vertical induced velocity along radius 
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Fig 34. Variation of Induced angle-of-attack along radius 
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\f!Uicai!on for a fixed wing 
( ’ 1.inj.’iiage program is used 


Dcfifiitian iif \*ariab!es used 

a : 

lift curve coefficient 

3\‘( ’ : 

chord 

Ar : 

area of the wing 

Croot : 

root chord 

Clip : 

tip chord 

CL : 

coefficient of lift 

ganinia : 

circulation of the vortices 

ganimal : 

circulation of lifting line 

gaiiimaln : 

new circulation value of lifting line 

N 

number of control points 

rep : 

locations of control points 

rho 

density 

rA : 

stations on the wing 

sweep : 

sweep angle 

S : 

span 

thetaO : 

angle of attack 

theta! 

twist rate 

u : 

speed of plane m/s 

V : 

induced velocity 


#incIude<stdio.h> 

#include<math.h> 

main() 

int ijjj,N,flag,flagl; 

float a,avC,Croot,Ctip,rho,S,thetal ,u,Ar,fl ,f2,D,L; 

float C[20],rA[20],rcp[20],thetaO[40],v[20],CL[40],gamma[20], 
gammal[20],gammaln[20],sweep; 

printf("Enter the value of lift curve slope\n"); 
scanf("%f',&a); 

printf("Enter the value of average chord\n"); 
scanf("%f’,&avC); 

printf("Enter the value of root chordNn"); 

scanf("%f',&Croot); 

printf("Enter the value of tip chord\n"); 

scanf(''%f',&Ctip); 

printf("Enter the number of control points'®"); 
scanf("%d",&N); 

printf("Enter the value of densityXn"); 
scanf("%f'.&rho); 



p!inTf("I'nti‘t the span^n"); 

Kcantt""';,,r.&S); 

print ft "{inter the rate of twist'in"); 
scanlt "•!i>f'.&thetal ); 
pnntft "Enter the plane speed in m/s \n''}; 
scanf("";if',&u); 
printff* Enter the sweep angle”); 
scanfC‘®-(>f'. &sweep); 
theta0ll}-.0.07; 
for(i 2;i< - 33;if +) 

! 

thetaOfi I -thetaO[i- 1 ]+0.0 1 ; 

! 

Ar S*avC; 
for(i“l;i<--N;j++) 

I 

rA[i]-S/N*(i-l); 

rcp{i]--rA[i]+S/(2*N); 

} 

rAfN+l]-S; 

for(i=I;i<=5;i-H-) 

{ 

C[i]=Croot-(Croot-Ctip)*2/S*(S/2-rcp[i]); 

I 

for(i=6;i<= 1 

{ 

C[i]=Croot+(Croot-Ctip)*2/S*(S/2-rcp[i]); 

> 

for(jj=lyj<=lyj++) 

for<i=I;i<=N;i++) 

{ 

gammal[i]=I; 

gainma[ 1 ]=-gammal[ 1 ]; 
gamma[N+ 1 ]=gammal[N]; 
for(i=2;i<=N ;i-+-+) 

{ 

gainma[i]=gaminal[i- 1 ]-gammal[i] ; 

} 

for(i=l ;i<=N;i-H-) 

{ 

v[i]=0; 

} 

for(i=l ;i<=N;i++) 

{ 

for(j=l y<=N+ 1 U++) 

v[i]=v[i]+gamma[j]/(4*22/7*cos(sweep)*(rA0]-rq)[i])); 

} 

} 

flag=l; 

while(flag=l) 

for(i=l ;i<=N/2;i-H-) 

{ 



B sqn(,fBfHO*f2); 

.i’*r|i]*B/2*(thetaO0i]-thetal* 

( S'.l-i cp[ i ) )5-atan{n /f2)+0.0262; 

} 

N;i++) 

I 

fl v{il; 
f2 

gamma!n[il=a*C[i]*f3/2*(thetaO[ij]+thetal* 

(S/2-rcp[i]))-atan(n /f2)+0.0262; 

f 

fiag!=0; 

for(i-l;i<=N;i++) 

I 

if^abs(gammaln[i]-gammal[i])>=1.0) 

{ 

flagl=l; 

} 

gammaI[i]=gammal[i]+0. 1 *(gammaln[i]-gaTnmal[i]); 

gamma[ 1 ]=-gammal[l ]; 
gamma[N+ 1 ]=gammal[N]; 
for(i=2;i<=N;i++) 

{ 

gamma[i]=gammal[i-l]-ganimal[i]; 

if(flagl=0) 

{ 

flag=0; 

break; 

} 

for(i=l;i<=N;i-H-) 

{ iff 

v[i]=0; ■ 

} 

for(i= 1 ;i<=N ;i++) 

{ 

for(j=l ;j<=N+l ;j++) 

{ 

v[i]=v[i]+gamma[j]/(4*22/7*cos(sweep)*(rA[j]-rcp[i])); 

} 

} 

} 

for(i=l ;i<=N;i++) 

L=L+rho*u*cos(sweep)’'‘gamtnal[i]*S/N; 

} 

CL|jj]= 2*L/(rho*u*u*S*avC); 

printf("thetaO = %f \t CL = %f\n", thetaO[jj],CL[jj]); 



I i«i a helicopter rotor blade using 
a picsciihed uake model (MAIXAB code) 

‘bi Ikfiiiition of Variables used 
a ' lift ciinc coefficient 
ft as : angle step size 
A : a constant 
•’* Ar : area of the blade 
: chord 

Cl) : coefficient of drag 
ti C 1 . : coefficient of lift 
% VQ : torque coefficient 
% C''r : Thrust Coefficient 
% dCT : differential Thrust Coefficient 
giimnia : circulation of the vortices 
% gaiiutia!: circulation of lifting line 
% gaimmn: new value of gamma 

gammat: total circulation of tip vortex 
% k I : Vortex Settling Rate Parameter 1 
% k2 : Vortex Settling Rate Parameter 2 
% !c : polynomial fitting lift-curve 
% lambda 1 : Wake Contraction Rate Parameter 
% lambda : total inflow ratio 
% lambdac: axial inflow ratio 
% lambdai: induced inflow ratio 
% mu : inverse lift to drag ratio 
% N : no. of blades 
% novs : no.of vortex stations 
% ntv : no. of tip vortices 
% omega : angular velocity of blade 
% p : polynomial fitting drag-polar 
% psib : angle between the blades 
% psiw : azimuth angle 
% R : radius of a blade 

% rl : 

% rA : stations on the blade 
% rb : non-dimensional radius 
% rc : root cut-off 
% rep : locations of control points 
% rho : density 
% rtif : root tip loss factor 
% Rul ; Roll up location of tip vortices 
% Rux,Ruy,Ruz : coordinates of the roll up 
% location 

% sigma : solidity ratio (constant chord) 

% sweep: tip sweep angle 

%T -.Thrust 

% Treq : Thrust required 

% theta : angle of attack 

% thetaO : collective pitch 

% thetal : twist rate (radians per unit length) 

% thetal R: total twist (in degrees) 

% V : axial velocity 
% vh : mean induced velocity in hover 
% vx,vy,vz : total induced velocity components 
% vxf,v;^,vzf :farwake induced velocity 



yiA/i inhodid vortex induced velocity 
‘'\i : lifting line induced velocity components 

\ : indiicetl velocity components of tip vortex 
''^1 xiAK/i : coofiiinateH of inboard vortices 
: 

xt,yi,/t : coordinates of tip vortex 
*^0 Bkiiie assumed to be initially along the Y-axis 

global a as A Ar C Cl. CD CQ CT fl f2 B gamma gammal 

global kl k2 Ic lambda lambdal lambdac 
global N iitv novs 
global omega p 

global psib R rA rep re rho rtif Rul Rux Ruy Ruz 

global sigma sweep 

global T theta 1 theta! R thetaO Treq V vh vx vy vz 

global vzt vzi vzf 


fidl - fopen(*bladeegl.nikY); 
a - fscanf(fidl,*%f,I); 
as-=fscanf(fidi;%f,l); 

C«fscanf(fidi;%f,}); 

N = fscanf(fidi;%f,l); 
omega - fscanf{fidU’%f,l); 
rho«fscanf{fidi;%fd); 

R = fscanf(fidi;%f,l); 
re « fscanf(fidl,’%f J); 
rtIf- fscanf(fidl,’%f,l); 

Rul«fscanf(fidi;%f,l); 

V = fscanf(fidi;%f,l); 
thetaO » fscanf(fidl ,'%f , 1 ); 
thetal=^fscanf(fidi;%f,l); 
sweep = fscanf(fidl,’%f ,1); 
status = fclose(fidl); 

A = 0.78; 

Ar=R*C; 
psib = 2*pi/N; 
sigma = N*C/(pi*R); 
lambdac = V/(omega*R); 
lambda = 0.; 
theeta = 0.; 

thetal R 1 80/pi*thetal *(1 -rc)*R; 

CL=0;CD=0;CQ=0;CT=0;fl=0;f2=0;B=0;kl=0;k2=0;lambdal=0 

Rux=0;Ruy=0;Ruz=0;T=0;vx==0;vy=0;vz=0; 


% Discretize the blade 
% vortex stations 

novs = round((rtlf-rc)/0.05) + 1 ; 
ntv = round((l-rtlf)/0.01); 
for i = 1 :novs 

rA(i) = (rc + (i-l)*(rtlf-rc)/(novs-l))*R; 
end 



fin I l:iilv 
tA|i*ii0vs|- rtlPR4 
end 

CiW!fcil points 
for I ‘ !;{iiovs*l ) 

fcp(i) - rA{i)Hrt!f-rc)/(2A(iiovs‘-l))*R; 
eiici 

for i ' " i :ntv 

rcp{i*^ novs*l) -rA(i+iiovs-l)+(I-rtIf)/(2.**'ntv)*R; 

end 

0 ^ is ?|{ >fg ♦♦ ♦ Ift ♦ J»t ]|t i»t J|t * sjc Ijj t ♦*♦**♦*♦♦♦**♦ * 

f 0 

% polynomial fitting lift-curve 

xl - [-0.2793 -0.2443 -0.2094 -0.1745 -0.1396 -0.1047... 

-0.0698 -0.0349 0.0 0.0349 0.0698 0.1047 0.1396... 

0.1745 0.2094 0.2443]; 

yl ■= [-1.6 -1.45 -1.3 -l.I -0.85 -0.625 -0.4 -0.2 0.0... 

0.225 0.4 0.625 0.8 0.925 1 .0 0.6]; 

Ic = spline(xl,yl); 

% polynomial fitting of drag-polar 
X = [-0.8 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0... 

0.1 0.2 0.3 0.4 0.5 0.6 0.8]; 

y = [0.018 0.0134 0.0121 0.01 12 0.0106 0.0102 0.0099 ... 

0.0098 0.01 0.0102 0.0108 0.01 14 0.0124 0.014 0.018]; 
p = spline(x,y); 


[gamma, gammal] = inigam(l); % initialisation of circulations 
[CT] = ThrCoe(theta0); 
vh = 1 .0; 

[ vx,vy,vz]=TotInd( 1 ); 

[gamma] = Compgam(l); 

[CT,CQ,alpha] = PerCal(l); 

% Plots of Performance calculations 

pIot(CT,CQ) 

figure(2) 

plot(CT/sigma,CQ/sigma) 

figure(3) 

plot(rcp/R,alpha) 

figure(4) 

plot(rcp/R,gammal) 

o/„*********K****»**i‘***fEnd of main program***************************** 

% function for initialization of the blade bound circulation(gammal) values 

function[gamma,gammal] = inigam(cc) 

global a as A Ar C CL CD CQ CT fl £2 f3 gamma gammal 

global kl k2 Ic lambda lambdal lambdac 

global N ntv novs 

global omega p 

global psib R rA rep rc iho rtlf Rul Rux Ruy Ruz 
global sigma 

global T thetal thetalR thetaO Treq V vh vx vy vz 



V/! v/i v/f 


C”T2 : (*T pcT imil length 
ticchic : deciciifig factor 
^’rtT2 : 1' per unit length 

thee!ii2 "(ihetaO‘t theta! '**(rcp'rc’*‘R)); 
ilecfac f!aiiiMac/2“i‘Sigma*a/I6)^2 + ... 

.sigma*a'S‘*‘(theeta2.’^rcp/R‘"lambdac); 
if aiiy(decfac<0.0) 
cm I :iKws4otv-l )- 0.003/R; 
else 

lambdail -- -(lambdac/2+sigma*a/I6)+ ... 

sqrt((laTnbdac/2+sigma*a/16)^2 + ... 
S!gma*a/8*(theeta2.’*'rcp/R-lambdac)); 

Iambda2 = lambdai2 + lambdac; 

CT2 « sigma*a/2*(theeta2,*(rcp/R).'^2"Iambda2.*rcp/R)*l/R; 

end 

T2 = CT2*rho*pim^2*(omega*Rr2; 

gaiiimal{! :novs+ntv-l ) = 0; 
gamma{ 1 :iiovs+ntv) =■ 0; 

gammal = T2./(rho*(omega*rcp)); 
for i ~ novsinovs+ntv-I 

gammal(i) = (rA(novs+ntv)-rcp(i))/(rA(novs+ntv)... 
-rcp(nDvs))*gamnnial(novs-l ); 
end 

gamma( 1 )«-gammaI( I ); 
gamma(no vs+ntv)=gammal(no vs+ntv- 1 ) ; 
for i =« 2:novs+“ntv-l 
gamma(i)«gammal(i*- 1 )--garnmal(i); 
end 

% ♦♦>»t*»»5»|.***>t«J|«*>*e^********************************’*'*************** 

% Function to estimate Thrust Coefficient 
function[CTl] = ThrCoe(thetaO) 

global a as A Ar C CL CD CQ CT f 1 f2 f3 gamma gammal 
global kl k2 Ic lambda lambdal lambdac 
global N ntv novs 
global omega p 

global psib R rA rep rc rho rtlf Rul Rux Ruy Ruz 
global sigma 

global T thetal thetal R thetaO Treq V vh vx vy vz 
global vzt vzi vzf 

% non-uniform inflow 
lambdac = V/(omega*R); 
dr =[rc:(L0-rc)/l 0:1.0]; 

CTl =0; 
fori = 1:10 

CTl = CTl +quadlCdCT;dr(i),dr(i+l)); 
end 



1 CTI •rlio*pi*R'*2^(onicga"^R)'^2; 


litiictHm for tliriist coefficient 
function f “* dCl^rb) 

global a as A Ar C CL CD CQ CT f 1 f2 f3 gamma gammal 
global kl k2 Ic lambda lambda I lambdac 
global N ntv novs 
global omega p 

global psib R rA rep re rho rtlf Ru! Rux Ruy Ruz 

global sigma 

global T theta! theta IR thetaO Treq V vh vx vy vz 

global vzt vzi vzf 

% decfac : deciding factor 
% rb : r/R 

theeta = {thetaO + theta 1 *rb*R); 
decfac ~ (lambdac/2+sigma*a/l6)^2 + ... 

sigma*a/8*(theeta.*rb-Iambdac); 
if any(decfac<0.0) 
f = 0,003’^tb./rb; 
else 

lambdai = -(iambdac/2+sigma*a/16)+ ... 

sqTt((lambdac/2+sigma*a/16)^2 + 
sigma*a/8*(theeta.*rb“lambdac)); 
lambda ~ lambdai + lambdac; 

f « sigma*a/2'*'(theeta.*(rb).'^2-lambda.*rb); 

end 

***’!'♦*♦**♦************** 

% Total Induced velocity calculation 
function[vx,vy,vz]==TotInd(cc) 

global a as A Ar C CL CD CQ CT f 1 f2 B gamma gammal 
global kl k2 Ic lambda lambdai lambdac 
global N ntv novs 
global omega p 

global psib R rA rep re rho rtlf Rul Rux Ruy Ruz 
global sigma 

global T thetal thetalR thetaO Treq V vh vx vy vz 
global vzt vzi vzf 

[vxt,vyt,vzt]=IndTip(l ); % tip-vortex contribution 

[vxi,vyi,vzi]=IndInbV(l); % inboard-sheet contribution 
[vxl,vyl,vzl]=IndLL(l); % other lifting-line contribution 
[vxf,vyf,vzq=IndFW{l); % far-wake contribution 

vx = vxt-l-vxi+vxl+vxf; 
vy = vyt-l-vyi+vyl+vyf; 
vz = vzt+vzi+vzl+vzf; 

o^*************************************************************** 



fo! calculating the induced velocity due to tip-vortex 

funiiiofi!vx!.v\t,\/t] IndTip{cc) 

jdofxil a as A Ai C’ Cl. CD CQ CT fl f2 fi gamma gammal 
f!la!>;i! kl k2 Ic lambda lambda! lambdac 
fjohal N fitv novs 
gk>hal cuBcga p 

global pstb R rA rep rc rho rtlf Rul Rux Ruy Ruz 

global sigma 

global 'F theta! thetal R tlietaO Treq V vh vx vy vz 
global vzt vzi vzf 

vxf{ 1 movS't'ntv-l )“0;vyt(l :novs+ntv-l)=0;vzt(l :novs+ntv-l)=0; 
% initialisation of induced velocities 


fork 

% Lifting Line equation 
% assumption -> reference blade on y axis 
% assumption -> Rul < psib 
% rotating coordinate system chosen 

gammat 0; 
for i - novs+1 :novs+ntv 
gammat = gammat + gamma(i); 
end 

for 1 « l:novs+ntv-I 
count = round{(Rul-0)/as); 

ang == [0+k*psib:(Rul-0)/count:Rul+k*psib]; % dividing the domain 
for jj « 1 mount 
for i = I :ntv 

vxt(l) « vxt(l) 4- gamma(i+novs)/(4*pi)*quadlCpx',ang0j),ang0j+l),[],[],i,Lk,l,l); 
%vyt(I) = v>t(I) + gamma(i+novs)/(4*pi)*quadl(’px’,ang(jj),ang(jj+l),[],[]d,2,k,l,l); 
^ vzt(l) = vzt(l) + gamma(i+novs)/(4*pi)*cos(sweep)*ang0j),ang0j+l),[],[],i,3,k,l,l); 
end 
end 

count = round((psib-Rul)/as); 

ang “ [Rul+k*psib:(psib-Rul)/count:psib+k*psib]; 

forjj = I ’.count 

vxt(l) ” vxt(l) 4- gammat/(4*pi)*quadlCpx’,ang(jj),ang(jj4*l),[],[],ntv,l,k,l,l); 

%vyt(l) = vyt(l) + gamTnat/(4*pi)*quadlCpxkang(jj),ang0j4-l),[],[],ntv,2,k,l,l); 
vzt(l) = vzt(l) 4- gammat/(4*pi)*cos(sweep)*ang(jj),ang0j4-l),[],[],ntv,3,k,Ll); 
end 

count = round((4*psib-psibyas); 

ang = [psib+k*psib:(4*psib-psib)/count:4*psib4-k*psib]; 
forjj = 1 mount 

vxt(l) = vxt(l) + gammat/(4*pi)*quadlCpxkang(jj),angOj4'l),[],[]>ntv,l,k,l,0); 

%vyt(l) = vyt(l) + gammat/(4*pi)*quadl(’px',ang(ij),ang(jj-M ),[],[], ntv,2,k,l,0); 
vzt(l) = vzt(l) 4“ gammat/(4*pi)*cos(sweep)*angGj),ang(jj+l),[],[],ntv,3,k,l,0); 

end 

end 



% iiitr|:ra!itl for tip vortex induced velocities 
functicni f p\(psiwj,j\kj,flag) 

global a as A Ar C CL CD CQ CT fl f2 O gamma gammal 
global kl k2 Ic lambda lambda I lambdac 
global N ntv no vs 
global omega p 

global psib R tA rep rc rho rtlf Rul Rux Ruy Ruz 

global sigma 

global T thetai theta! R thetaO Treq V vh vx vy vz 
global vzt vzi vzf 

% i : tip vortex no. 

%j : component (x«!»y=2,z=3) 

% k : blade no. 

% 1 : control pt. no. 

% flag : just a no, 

(Rux, Ruy, Ruz] = TipVor(Rul+k*psib,k); 


if is=*«ntv 

[xt(i.;),yt(i,:), 2 t(i,;)] = TipVor{psiw,k); 

else 

[xt(:,:),yt(:,:), 2 t{:,:)] = OthTV(psiw,k,i); 
end 

% Lifting Line equation 
% assumption -> reference blade on y axis 
% assumption -> Rul < psib 
% rotating coordinate system chosen 

if (i~===mtv)&(psiv^KRtil+k^ 
dxt(i,:) = (-l/Rul)*rA(i+novs)*(-sin(k*psib))+Rux/Rul; 
dyt(i,:) = (-l/Rul)*rA(i+novs)*cos(k*psib) + Ruy/Rul; 
dzt(i,:) = Ruz/Rul; 
else 

dxt(i,:)= •-R*(-lambdal)*(l-A)*exp(-lambdal *(psiw-k*psib))... 
.*sm(psiw) - R*(A+(l-A)*exp(-lambdal*(psiw-k*psib)))... 
.*cos(psiw) ; 

dyt(i,:)= R*(-lambdal )*(1 -A)*exp(-lambdal *(psiw-k*psib))... 

.*cos(psiw) - R*(A+(l-A)*exp(-lambdaP(psiw-k*psib)))... 
.*sm(psiw) ; 
dzt(i,:>=<R*kl); 
end 

if flag == 0 
dzt(i,:)=-(R*k2); 
end 


rl(i,:) = sqrt((-rcp(l)*sin(0)-xt(i,:)).^2+... 

(rcp(l)*cos(0)~yt(i,:)).-^2+Czt(i,:)).^2); 



% to prevent singularity 


len l-'iitilhfrlCj,;!), 
for lii I:Ie!i 
if rUijiiK' 0.005*R 
!l|ijnl4f005*R; 
end 
eiiil 

Hvritch J 

case I , 

f -Cd>t(i,:).*zt(i,:)+dzt{i,:),*(rcp(l)*cos(0)-yt(i,:)))./(ri(i,:).^3); 

case 2, 

f ' (dzt(ia).‘*‘(“rcp(I)*sin(0)-xt(i,:))+dxt(i,:)-*2t(i,:))./(rl(i,:),'^3); 

case 3* 

f (dxt(i,:).'*‘{rcp(l)*cos(0)-yt(i,:))-dyt(i, . 

(*rcp(l)*sin(0)-xt{i,:)))./(rl{i,:)."^3); 

end 

0 / 1 ^ s^E ^ it if< s|( ^ 4 3|e 9|c ^ 4c 4c 4c 9tc 3|c ^ ^ ^ ^ ^ ^ ^ ^ >i<: 4c ’{c ^ ^ ^ ^ ^ sic ^ ^ He ’it ^ ^ ^ ^ 


% Function defining tip vortex geometry 
% Parametric equations 

fanction[xt5,yt5,zt5] = TipVor(psiw,k) 

global a as A Ar C CL CD CQ CT fl f2 O gamma gammal 
global k! k2 Ic lambda lambdal lambdac 
global N ntv novs 
global omega p 

global psib R rA rep re rho rtlf Rul Rux Ruy Ruz 
global sigma 

global T theta 1 theta 1 R thetaO Treq V vh vx vy vz 
global vzt vzi vzf 

thetalR ~ thetal*(l-rc)’*'R*180/pi; 

A « 0.78; 

B--0.000729*thetalR; 

Cl »2.3 + 0.206*thetalR; 
m “ 1.0 - 0.25*exp(0.040*thetalR); 
n = 0.5-0.01 72*thetalR; 

kl = (B + Cl*(CT/NAtirm )*(L35-0.3/exp(V/vh)); 

CTO - N^n*(-B/Cir(l/m); 

k2 = -(abs(CT-CT0)r0.5*(1.35-0,3/exp(V/vh)); 

lambdal = 4.0*(CT)^0.5*(-0.35+0.47*(1.7+V/vh)^2); 

xt5 = -R*(A+(l-A)*exp(-lambdal*(psiw-k*psib)))-*sin(psiw); 
yt5 = R*(A-J-(l-A)*exp(-lambdal*(psiw-k*psib))).*cos(psiw); 

if (psiw>~(0+k*psib)&psiw<=(psib+k*psib)) 
zt5 == -(R*kl*(psiw-k*psib)); 
elseif (psiw>~(psib+k*psib)&psiw<=(4*psib4-k*psib)) 
zt5 == -(R*(kl*psib+k2*(psiw-k*psib-psib))); 
else 
zt5 = 0; 
end 



t(pm\ '' 4*psil'i*k'^psih 
\t5 Cl; \i5 0; zt5 ' 0; 

eiiii 


i> t' ft ♦♦♦in# ♦♦♦»!< ♦♦«!«♦♦ + **♦**»}<♦♦*% ****>!«* 


% Function defining tip vortex geometry 
^Uf Parametric equations 

fiiiictiO!ifxt5,>l5,zt5] - TipVor(psiw,k) 

global a as A Ar C CL CD CQ CT f 1 12 f3 gamma gammal 

global kl k2 Ic lambda lambdal lambdac 
global N ntv novs 
global omega p 

global psib R rA rep rc rho rt!f Rul Rux Ruy Ruz 

global sigma 

global T theta! theta! R thetaO Treq V vh vx vy vz 

global V2t vzi vzf 

theta! R = theta! *(l-rc)*R* 1 80/pi; 

A = 0.78; 

B«-0.000729*thetalR; 

Cl «-23 + 0.206*thetaIR; 
m ^ 1.0 - 0.25*exp(0.040*thetalR); 
n = 0.5-0.0172*thetaIR; 

kl «« (B + C1*(CT/N^nrm )*(1.35-0.3/exp(V/vh)); 
CTO^N^rf^C-B/CirCl/m); 
k2 « «(abs(CT-CT0)r0.5^(1.35-0.3/exp(V/vh)); 
lambdal « 4.0‘'(CTr0.5*(-035+0.47*(l .7+V/vhr2); 

xt5 « -R*(A+(l-A)’*‘exp(-lambdal’‘'(psiw-k*psib))).*sm(psiw); 
ytS - R*(A+'(l"A)’*'exp(-lambdal*(psiw-k*psib))).*cos(p$iw); 

if (psiw>=(0+k*psib)&psiw<=(psib+k*psib)) 
zt5 = -(R*kl*(psiw-k*psib)); 
elseif (psiw>"(psib+k*psib)&psiw<=(4*psib+k*psib)) 
zt5 = -(R*(kl*psib+k2*(psiw-k*psib-psib))); 
else 
zt5 = 0; 
end 


if psiw > 4*psib+k*psib 
xt5 = 0; ytS = 0; zt5 == 0; 
end 




% Equations of other tip vortices 

function[xt,yt,zt] = OthTV(psiw,k,i) 

global a as A Ar C CL CD CQ CT fl fZ f3 gamma gammal 
global kl k2 Ic lambda lambdal lambdac 
global N ntv novs 
global omega p 



|»!ohai psih R lA rep re rho illf Riil Rux Ruy Ruz 
pjp|i;il sipma 

1“ tlicial theta I R thetaO 1’req V vh vx vy vz 
pjiibal \/l v/f 

t |psi\v4*psib)'Riil; 

ifpsiw< (Rul 4 k*psib) 

C I «t)*rA(i“tnovs)*(-si!i(k*psib)) + t*Rux; 

) { 1 ••i)*rA(i4“novs}*cos(k*psib) + t'*'Ruy; 

zRiel ' t*Riiz; 

elseif {psi\v>{Ru!4'k^psib)&psiw<-(psib+k*psib)) 

"* »R^(A+(l-A)*exp(-lambdal*(psiw-k*... 
psib))).*sin(psiw); 

yt(i,:) - R*(A4(l-A)’^exp(-lambdaI*(psiw-k*... 

psib))).*cos(psiw); 

2 t(i,:) =" *(psiw'-k*psib)); 

elseif (psiw>(psib+k*psib)&psiw<=(4*psib+k*psib)) 
-R*(A+(I-A)*exp(“lambdal*(psiw-... 
k*psib))).*sin(psiw); 

= R’^(A+(l-A)*exp(“Iambdal*(psiw-... 
k*psib))).*cos(psiw); 

2t(i,:) = + k2*(psiw-k*psib>psib))); 

else 

- 0 ; 

yt(i,:)-0; 

= 0 ; 

end 

% Induced velocities due to Inboard Vortices 

function[vxi,vyi,vzi] - IndlnbV(cc) 

global a as A Ar C CL CD CQ CT fl f2 f3 gamma gammal 
global kl k2 Ic lambda lambda 1 lambdac 
global N ntv novs 
global omega p 

global psib R rA rep re rho rtlf Rul Rux Ruy Ruz 
global sigma 

global T thetal thetalR thetaO Treq V vh vx vy vz 
global vzt vzi vzf 

vxi(l :novs+ntv-l)=0;vyi(l :novs+ntv-l)=0;vzi(l :novs+ntv-l )=0; 
if psib>pi/2 

angl = pi/2;ang2 = psib; 
else 

angl = psib;ang2 == pi/2; 
end 

fork = 0:N-1 
countl = round(angl/as); 
anga = [04“k*psib:angl/countl:angl+k*psib]; 
count2 = round((ang2-angl)/as); 



if fiingl iin*!2) 

1 :c(>ii!it2) » Cl; 

eke 

a«|»h • jiiiigl ‘ i^’^psib'fang2-afigl)/coiint2:ang2‘l“k*psib]; 

end 

coiml? ' !ouiKl|'f4'*‘psib-ang2)'^as); 

aiigc |ang24'k*psib:{4*psib-ang2)/coi2nt3:4*psib+k*psib]; 

fiir ! ■ I 

forjj Ixoimtl 
for i 1 :iiovs 

vxi(!) vxi(!) + gamma(iy(4*pi)*quadl('qx>nga(jj),anga(jj+l), [],[], i,l,k,l,0); 
o«vyi(l) = vyi{I) + gamma(i)/(4*pi)*quadl(’qx’,anga(jj),anga(jj+l),[],[],i,2,k,l,0); 
V2i(l) ™ V2i(l} + gamma(i)/(4*pi)*quadl(’qx',anga(jj),anga0j+l),[],[i,i,3,k,l,0); 

end 

end 

forjj = l:count2 
for i = I :novs 

vxi(l) = vxi(I) + gamma(i)/(4*pi)*quadl('qx',angb(jj),angb(ij+l),[],[],i,l,k,l,I); 
%vyi(l) = vyi(l) + gamina(i)/(4*pi)*quadl('qx',angb(jj),angb(jj+l), [],[], i, 2, k, 1,1); 
vzi(l) = vzi(l) + gamma{i)/(4*pi)*quadl('qx',angb(jij),angb0j+l), [],[], i, 3, k,U); 

end 

end 

forjj = l:count3 
for i “ 1 :novs 

vxi(l) = vxi(l) + gamma(i)/(4*pi)*quadl('qx'.angc(jj).angc0j+l),[],[],i,l,k,l,2); 
%vyi(l) = vyi(l) + gamma(i)/(4*pi)*quadl('qx',angc(iij),angc(ij+l),[],[],i,2,k,l,2); 
vzi(l) = vzi(l) + gamma(i)/(4*pi)*quadlCqx',angc(jj),angc(jj+l),[],[],i,3,k,l,2); 
end 
end 
end 
end 


3ft 4t ;%i 3)1 :<( Jj|t 9fc :ie 3|e ]fc 3)c 3fe jfc 9(1 % ;:)e IK 4c * ♦ * * % 3ic 3)c :4c 9fc He * 3|c 4: * 3ie 4c * :>ic ^ 
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% integrand for inboard vortex induced velocities 
function f = qx(psi,i j,k,I,flag) 

global a as A Ar C CL CD CQ CT fl f2 f3 gamma gammal 
global ki k2 Ic lambda lambdal lambdac 
global N ntv novs 
global omega p 

global psib R rA rep rc rho rtlf Rul Rux Ruy Ruz 
global sigma 

global T thetal thetalR tbetaO Treq V vh vx vy vz 
global vzt vzi vzf 

% i : vortex no. 

%j : component (x=l,y=2,z=3) 

% k : blade no. 

% I : control pt. no. 

% flag : just a no. 



Uh llkialEi2K*|(l.45*the^alR-^I8))*sqrt(CT^^^ 
kl 1 -2,r sqiKl ' r, ?r( 1 .354).3/exp(V/vh)); 

k2 1 ^*2 ?*sqiliC"*lV2l*( 1 .35-0.3/exp(V/vh)); 

switch 
case fh 

7HI/J “(kl I •(psi-k*psib)*rA(!)); 
ci/!(ial -(kl i*rA(i)): 


case 1, 

if psib pi/2 

" “(k20*(psi-k*psib-pi/2)*R+(kl 1*... 
(psi“k*psib}"k20’^(psi“k**'psib... 
-pi/2)frA(i)); 

- .(k20*R+(kl i-k20)*rA(i)); 

else 

zi{i,:) “ -(kl 1 *psib4-k2I *(psi-k*psib-psib)}... 

^rA(i); 

d7j{i,:)«-{k2PrA(i)); 
end 
case 2 

2i(ic-) = "(R’^k20’^(psi“k*psib-pi/2)+... 

(kl 1 ’*‘psib+k21 *(psi-k*psib-psib)-... 
k20*(psi-k*psib-pi/2))*rA(i)); 
dziCi,:) - -(R^k20+(k21-k20)*rA(i)); 
end 


if zi(i,:) *== 0.0 
psiw = 0.0; 

elseif (zi(i,:)<-(kl ♦psib*R))&(zi(i,:)'-=0.0) 
psiw » -zi(!,:)/(R*kl)+k*psib; 
else 

psiw = (-zi(i,:)/R-kl *psib)/k2 + psib + k*psib; 
end 

if psiw > 8*pi/N+k*psib 
psiw = 8*pi/N+k*psib; 
end 

rt = R*(A+(l-A)*exp(-lambdal*(psiw-k*psib))); 
rad == rA(i)*rt/R; 
xi(i,:) = -rad.*sin(psi); 
yi(i,:) = rad.*cos(psi); 

dxi(i,:) = -rad.*cos(psi); 
dyi(i,:) = "rad.*sin(psi); 

r2(i,:) = sqrt((-rcp(l)*sin(0)-xi(i,:)).^2+(rcp(l)*cos(0)-yi(i,:)).^2.... 
+(-zi(i,:)i^2); 

len = length(r2(i,:)); 
foriii = l:len 
ifr2(i,iii)< 0.005*R 
r2(i,iii)= 0.005*R; 
end 
end 


% to prevent singularity 



switch j 
case 1 , 

case 2, 


f (tii?i(i.:).*(-rcp{I)*sin(0)-xi(i,;))+dxi(i,:).*zi(i,:))./(r2(i,:).''3); 

ease 3* 


f <d';i(i.;).*(rcp(l)*cos(0)-yi(i,:))-dyi(i,:).*(-rcp(l)*sin(0)-xi(i,;)))./(r2(i,:).''3); 

eiid 




Induced velocities due to other Lifting Lines 


tiinct!Oiifvx!,vyhvzI]=IndLL(cc) 

global a as A Ar C CL CD CQ CT f 1 f2 f3 gamma gammal 
global k! k2 Ic lambda lambda 1 iambdac 
global N ntv novs 
global omega p 

global psib R rA rep rc rho rtlf Rul Rux Ruy Ruz 
global sigma 

global T theta 1 theta i R thetaO Treq V vh vx vy vz 
global V2t vzi vzf 


vxl(l :novs+ntv-l )-0;vyl(l :novs+ntv-l)=0;vzl(l :novs+ntv-l)=0; 

fork* 1:N-1 
for I * l:novs+ntv-l 
fori * l:novs+ntV"i 

V2l(l) «= vzl(l)+gammal(i)/(4*pi)*quadl('sx',rA(i),rA(i+l ),[],[].k,l); 
end 
end 
end 

% note : X &. y components are zero 

4c 4c a|c 4c 9{c 4c 4c >ie * sic 4c 4e 4c 41 ; 4e 4e 4c 4 k 4t 4e 4e 4e 4: 9|c 4c 4c 4e 4e 4c % 4c ^ % He 4: ’ic 4: >)( 4: 4c 4c 4c % 4: 4: 4e 4c 4: 4; 4: 4: 4t 4: 4e 4: 4e 4; % 


% integrand for lifting line induced velocities 
function f = sx(m,k,I) 

global a as A Ar C CL CD CQ CT fl f2 O gamma gammal 
global kl k2 Ic lambda lambdal Iambdac 
global N ntv novs 
global omega p 

global psib R rA rep rc rho rtlf Rul Rux Ruy Ruz 
global sigma 

global T thetal thetal R thetaO Treq V vh vx vy vz 
global vzt vzi vzf 


% k : blade no. 

% 1 : no. of control pt. 



rl - sqr1((-rcp(l)*sin(0)+m*sin(04k*'psib)).''2 + ... 
(^cp(l)*cas(0)-m*cos(0^k♦pslb)).''2); 


dix - - -sin(0 * k*psib); 
dly ■ cosCn ' k'*'psib); 
dIz-0; 

“/o X & y components of ind vel = 0 

f - Cdlx.*(rq>(l)*cos(0)-m*cos(0+k*psib))-... 

dly.*(-rcp{l)*sin(0) !-m*sin(0+k’''psib)))./(rl .'^3); 

fi / ♦ III Ifk 41 ]fi ]|( iE ^ i|r 4t sfe % sfc 94( % )ie ifc 4c ^ lie )|C * 4c 4; % 9)c He * * * * * * * * 9ie He % ai; 


% Induced velocity due to Far Wake 

fuiictioii[vxf,vyf,vzf]==IndFW{cc) 

global a as A Ar C CL CD CQ CT fl O O gamma gamma! 
global kl k2 Ic lambda lambda! lambdac 
global N ntv novs 
global omega p 

global psib R rA rep rc rho rtlf Rul Rux Ruy Ruz 
global sigma 

global T theta! theta I R thetaO Treq V vh vx vy vz 
global vzt vzi vzf 

vxf(l:iiovs+ntv-l) = O;vyf(l:novs+ntv-l) = O;vzf(i:novs+ntv-l) = 0; 
gammat = 0; 
for j « novs+1 :novs+ntv 
gammat = gammat + gamma(i); 
end 

count = round(6.2832/a$); 
ang = [0:6.2832/count:6.2832]; 

forjj =1 .'count 
for I = 1 :novs+ntv-l 

vxf(l) == vxf(l) + 8*gammat/(4*pi)*quadl(Vx’,ang(jj),ang(jj-Hl ),[],[], Id); 
%vyf(l) = vyfO) + gammat/pi*quadl(’rx\ang(jj),angOj+l ),[],[], 2d); 
vzf(l) = vzf(l) + 8*gammat/(4*pi)*quadlCrx',ang(jj),ang(jj+l),[],[]»^d); 
end 
end 


% integrand for far wake induced velocities 
function f = rx(phi jd) 

global a as A Ar C CL CD CQ CT f 1 f2 f3 gamma gammal 
global kl k2 Ic lambda lambda 1 lambdac 
global N ntv novs 
global omega p 

global psib R rA rep rc riio rtlf Rul Rux Ruy Ruz 
global sigma 

global T theta 1 thetalR thetaO Treq V vh vx vy vz 



global vzt vzi vzf 


: component {X“ I,}^*2,2~3) 

1 : ctHitrol pt iK>, 

fxth.vlhjthl *' TipVor(4*psib,0); 


11- 4.0nCT)"0.5; 

ni “R*( A+“( 1 -A)*exp(-Il *8’*'pi/N)); % radius of the tipvortex in hover 

12- 4.0*(CT)"^t).5*(-0.35+0,467*(I.7+V/vh).^2); 

rr2“R*(A ^(l-A)*exp(-!2*8*pi/N)); % radius in axial flight 

xf == -iT2/iTl*R*sin(phi4-0); 
yf = trlirrl *R*cos(plii+0); 
zf = 2 th; 


rlx(l,:) = -rcp(l)*sin(0)-xf; 
r!y(l:) = rcp(l)^cos(0}-yf; 
r! 2 (l,:) = -zf; 

- sqrt(r!x(h:),^2 + rl><h:).^2 + riz(l,:).^2); 

dlx = -rr2/rrl ’^R*cos(phi+0); 
dly = -rr2/rrl *R*sin(phi+0); 
dlz «= 0; 
switch j 
case I, 

f=(dly.’^rlz(l/0"dl2.*rly(}g))./(rl(l,:).^3); 
case 2 

f = (dlz.*rlx(l,:)-<^l>^-’*'t'l2(h:))./(rl(l,:).'^3); 
case 3, 

f = (dlx.*rly(lg)-dly.*rlx(l,:))./(rl(l,:).^3); 
end 


♦**♦*♦%*♦**♦**♦♦**♦*♦*♦♦****♦*♦ ************** 


% Computation of gamma values 

function [gamma] = Compgam(cc) 

global a as A Ar C CL CD CQ CT fl f2 B gamma gammal 
global kl k2 Ic lambda lambda 1 lambdac 
global N ntv novs 
global omega p 

global psib R rA rep rc rho rtlf Rul Rux Ruy Ruz 
global sigma 

global T thetal thetal R thetaO Treq V vh Vx vy vz 
global vzt vzi vzf 

flagg = 1 ;countt=0; 
while flagg = 1 

for i = 1 :novs+ntv-l 
u(i) = omega*rcp(i); 
fl = V+vz(i); 
f2--(-u(i)+vx(i)); 
f3 = sqrt(fl.^2 + f2.^2); 



r.ininialnfi) an;*O/2.*(theta0+tlietal ♦(rcp(i)-rc*R>atan(fl^^ 

end 

tiaggi 0; 

for i l:novs4!itv-l 

if ahs(ganin]aln(i)»gammal{i))>(0.1 *gaimmaln(i)) 

llaggl - 1; 
end 

ganiitial(i) gamma}(i)+OJ*(gamma!n(i)--gamniaI(i)); 
end 

ganima( ! ) --gammaK I ); 
gamma(no vs+ntv)=^gamnial(novs+ntv- 1 ); 
for i = 2:iiovs+ritv-l 
gaiTima(i)=garnmal(i»l)"ganirnal(i); 
end 


if flaggl == 0 

flagg = 0; break; 
end 

[vx,vy, vz]=TotInd( 1 ); 
end 

3^*)|t** sit Ht*#***’fr**%>|c**** ♦*******************’(:*** 


% Performance calculations 

function[CT,CQ, alpha] = PerCal(cc) 

global a as A Ar C CL CD CQ CT f 1 f2 O gamma gammal 
global kl k2 Ic lambda lambdal lambdac 
global N ntv novs 
global omega p 

global psib R rA rep rc rho rtlf Rul Rux Ruy Ruz 
global sigma 

global T thetal thetalR thetaO Treq V vh vx vy vz 
global vzt vzi vzf 

% Calculate angle of attack 

alpha(l :novs+ntv-l) = 0; 

fori = l:novs+ntV“l 
u(i) = omega*rcp(i); 
n - V+vz(i); 
f2 =-(-u(i)+vx(i)); 
f3 = sqrt(fl.^2 + f2.^2); 
alpha(i) == 2*gammal(i)/(a*C*0); 

CL(i) = ppval(lc,alpha(i)); 

CD(i) = ppval(p,CL(i)); 
mu(i) = CD(i)/CL(i); 
drcp(i) = rA(i+l)-rA(i); 
gammab(i) = gammal(i)/(omega*R""2); 
end 

suml = 0; sum2 = 0; 
fori = l:novs+ntv-l 



Mini! * t41)fii*|V'‘\7(iV)Toniega*R)),.. 

» t f > * Ilf 1 1 1 i wiega* R'^2)*drcp(i)/R; 

Mini? *-11111? * \'\fi,r(amcga’*'R))) . 

n,l {*inkManr2)^rcp(i)/R^drcp(i)/R; 

elac 

Mini! Mini! ^ iicpliV'R“Vx(i)/(omega*R)»rnu(i)*(V+vz(i)).,. 

u.nu Vi^aannnaKi')*drcp(i)('R; 
siim? Mim? U\''’v/fU)'^>nicga*R)t*mu(i)*(rcp(i)/R-vx(i)... 
u>n!r]M*R in*panimab(i)’*'rcp(iVR*drcp(i)/R; 

end 

aid 

(?T N ; 

C*0 N'|ii*siim2; 



^ 4 



